function t = depcorrt(X,Y,Z,nZ)
% depcorrt  - dependent correlation comparison t statistic
%
% FORMAT:       t = depcorr(X, Y, Z [, nZ])
%
% Input fields:
%
%       X           vector (or matrix/array) X
%       Y           vector (or matrix/array) Y
%       Z           vector (or matrix/array) Z
%       nZ          optional vector (or matrix/array) nZ
%
% Output fields:
%
%       t           t-statistic for difference of correlations between
%                   rXZ and rYZ, d.f. = N - 3
% or
%                   z-statistic for difference of correlations between
%                   rXY and rZnZ
% 
%
% Note: the correlation is computed along the LAST non-singleton dimension
%
% See also: http://ssc.utexas.edu/consulting/answers/general/gen28.html
%     and:  http://luna.cas.usf.edu/~mbrannic/files/regression/corr1.html

% Version:  v0.7f
% Build:    8110521
% Date:     Nov-05 2008, 9:00 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 3 || ...
   ~isnumeric(X) || ...
   ~isnumeric(Y) || ...
   ~isnumeric(Z) || ...
    isempty(X) || ...
   ~isequal(size(X), size(Y)) || ...
   ~isequal(size(X), size(Z))
    error( ...
        'BVQXtools:BadArgument', ...
        'Invalid argument (combination) given.' ...
    );
end
if nargin > 3 && ...
    strcmp(class(nZ), class(X)) && ...
    isequal(size(nZ), size(X))
    v4 = true;
else
    v4 = false;
end

% translate?
sz = size(X);
if numel(sz) == 2 && ...
    sz(2) == 1
    ts = [1, 1];
    X = X';
    Y = Y';
    Z = Z';
    if v4
        nZ = nZ';
    end
    n = sz(1);
else
    ts = sz(1:end-1);
    n = sz(end);
    dm = prod(ts);
    X = reshape(X, [dm, n]);
    Y = reshape(Y, [dm, n]);
    Z = reshape(Z, [dm, n]);
    if v4
        nZ = reshape(nZ, [dm, n]);
    end
end

% three variables
if ~v4

    % compute the three correlations
    [rxy{1:2}] = cov_nd(X, Y);
    rxy = rxy{2};
    [rxz{1:2}] = cov_nd(X, Z);
    rxz = rxz{2};
    [ryz{1:2}] = cov_nd(Y, Z);
    ryz = ryz{2};

    % compute determinant
    % |R| =  (1 - rxy ^2 - rxz^2 - ryz^2 + (2*rxy*rxv*rvy)), the determinant of
    % the correlation matrix for X, Y, and V. 
    R = 1 + (2 .* rxy .* rxz .* ryz) - (rxy .* rxy + rxz .* rxz + ryz .* ryz);

    % compute dependent t-score
    t = (rxz - ryz) .* sqrt((n-1) .* (1 + rxy)) ./ ...
        (sqrt(2 .* ((n-1) / (n-3)) .* R + ((rxz + ryz) ./2 ) .^2 .* (1 - rxy) .^ 3));

% 4 variables
else
    
    % compute the required correlations
    [rxy{1:2}] = cov_nd(X, Y);
    rxy = rxy{2};
    [rxz{1:2}] = cov_nd(X, Z);
    rxz = rxz{2};
    [ryz{1:2}] = cov_nd(Y, Z);
    ryz = ryz{2};
    [rxn{1:2}] = cov_nd(X, nZ);
    rxn = rxn{2};
    [ryn{1:2}] = cov_nd(Y, nZ);
    ryn = ryn{2};
    [rzn{1:2}] = cov_nd(Z, nZ);
    rzn = rzn{2};
    rxyzn = (rxy + rzn) ./ 2;
    hxyzn = .5 .* ( ...
        ((rxz - rxy .* ryz) .* (ryn - ryz .* rzn)) + ...
        ((rxn - rxz .* rzn) .* (ryz - rxy .* rxz)) + ...
        ((rxz - rxn .* rzn) .* (ryn - rxy .* rxn)) + ...
        ((rxn - rxy .* ryn) .* (ryz - ryn .* rzn)));
    sxyzn = hxyzn ./ ((1 - rxyzn .* rxyzn) .^ 2);
    rxy(isnan(rxy)) = 0;
    rzn(isnan(rzn)) = 0;
    zxyzn = fisherr2z(rxy) - fisherr2z(rzn);
    t = sqrt(n - 3) .* zxyzn ./ sqrt(2 - 2 .* sxyzn);
    
end

% reshape to result
if numel(ts) < 2
    ts(2) = 1;
end
t = reshape(t, ts);
